Phase ordering and roughening on growing films 
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We study the interplay between surface roughening and phase separation during the growth of 
binary films. Already in 1+1 dimensions, we find a variety of different scaling behaviors, depending 
on how the two phenomena are coupled. In the most interesting case, related to the advection of a 
passive scalar in a velocity field, nontrivial scaling exponents are obtained in simulations. 
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Thin solid films are grown for a variety of technolog- 
ical applications, using molecular beam epitaxy (MBE) 
or vapor deposition. In order to create materials with 
specific electronic, optical, or mechanical properties, of- 
ten more than one type of particle is deposited. When 
the particle mobility in the bulk is small, surface configu- 
rations become frozen in the bulk, leading to anisotropic 
structures that reflect the growth history, and are dif- 
ferent from bulk equilibrium phases [Q. Characterizing 
structures generated during composite film growth is not 
only of technological importance, but represents also an 
interesting and challenging problem in statistical physics. 

In this paper, we examine the growth of binary films 
through vapor deposition, and study some of the rich 
phenomena resulting from the interplay of phase separa- 
tion and surface roughening. Simple models for layer by 
layer growth assume either that the probability that an 
incoming atom sticks to a given surface site depends on 
the state of the neighboring sites in the layer below [2], 
or that the top layer is fully thermally equilibrated [3|. 
Assuming that the bulk mobility is zero, once a site is oc- 
cupied, its state does not change any more. If the growth 
rules are invariant under the exchange of the two particle 
types, the phase separation is in the universality class of 
an equilibrium Ising model. Correlations perpendicular 
to the growth direction are characterized by the critical 
exponent v of the Ising model, and those parallel to the 
growth direction by the exponent ^£ m , with z m being the 
dynamical critical exponent of the Ising model. 

However, the layer by layer growth mode underlying 
these simple models is unstable, and the growing surface 
becomes rough. In many cases the fluctuations in the 
height /i(x, t), at position x and time t are self-affine, 
with correlations 
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where x is the roughness exponent, and Zh is a dynamical 
scaling exponent. A computer model with local sticking 
probabilities that allows for a rough surface was intro- 
duced in [Q. In 1+1 dimensions, the authors find phase 
separation into domains (with sizes consistent with the 
Ising model), and a very rough surface profile with sharp 
minima at the domain boundaries. We may ask the fol- 



lowing questions: (1) Are the roughness exponents dif- 
ferent at the phase transition point? (2) Are the critical 
exponents modified on a rough surface? We shall demon- 
strate that the coupling of roughening and phase sepa- 
ration leads to a rich phase diagram, and to nontrivial 
critical exponents already in 1+1 dimensions. 

To characterize phase separation, we introduce an or- 
der parameter m(x, i), which is the difference in the den- 
sities of the two particle types at the surface at position 
x and time t. The interplay between the fluctuations in 
m, and the height h is captured phenomenologically by 
the coupled Langevin equations, 
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d t m = K(V 2 m + rm — um 3 ) + aVh ■ Vm + bmV 2 h 
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Here, we have included the lowest order (potentially rel- 
evant) terms allowed by the symmetry m — > —m. Equa- 
tion (|2|) is the Kardar-Parisi-Zhang (KPZ) equation || 
for surface growth, plus a coupling to the order pa- 
rameter. Equation (||) is the time dependent Landau- 
Ginzburg equation for a (non-conserved) Ising model, 
with three different couplings to the height fluctuations. 
The Gaussian, delta-correlated noise terms, (h and £ m , 
mimic the effects of faster degrees of freedom. A differ- 
ent set of equations was proposed by Leonard and Desai 
H for phase separation during MBE. Their equations 
reflect the MBE conditions of random particle deposi- 
tion (in contrast to sticking probabilities that depend on 
the local environment), and a conserved order parameter 
which evolves by surface diffusion. They do not include 
the KPZ nonlinear ity. Computer simulations of corre- 
sponding 1+1 dimensional systems are presented in ||f7|. 

Dimensional analysis indicates that the couplings ap- 
pearing in Eqs. (^|3|) are relevant, and may lead to new 
universality classes. We shall leave the renormalization 
group analysis of these equations to a more technical pa- 
per, and focus here instead on computer simulations in 
1+1 dimensions. The quantities evaluated in the com- 
puter simulations are the height correlation function in 
Eq. (FA), and the order parameter correlation functions 
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perpendicular and parallel to the growth direction. Al- 
lowing for the possibility of different dynamic exponents, 
z m and Zh, for the order parameter and the height vari- 
ables, we fit to the scaling forms 

G$(x-x f ) = {m(x,t)m(x',t)) 

= \x-x'r 1 gi(\x-x'\/0 
G m \t-t') = (m(x,t)m(x,t')) 

= \t-t f \^ z -gi(\t-t f \/^). (4) 

Our simulations were done using a "brick wall" re- 
stricted solid-on-solid model (see Fig. [l]). Starting from 
a flat surface, particles are added such that no overhangs 
are formed, and with the center of each particle atop the 
edge of two particles in the layer below. We use two types 
of particles, A and B (black and grey in the figures). The 
probability for adding a particle to a given surface site, 
and the rule for choosing its color, depend on the lo- 
cal neighborhood. When A particles are more likely to 
be added to A dominated regions, and vice versa, the 
particles tend to phase separate and form domains. In 
this case, the order parameter correlation length £ is of 
the order of the average domain width. By changing the 
growth rules, it is possible to study cases in which some 
(or all) of the couplings a, 6, c, and a vanish, and thus 
to gain a complete picture of the different ways in which 
the height and the order parameter influence each other. 




FIG. 1. The "brick wall" model used in the simulations. 

The decoupled case, a = a = 6 = c = 0, is imple- 
mented using the following updating rules: A surface site 
is chosen at random, and a particle is added if it does not 
generate overhangs. Its color is then chosen depending on 
the colors of its two neighbors in the layer below. If both 
neighbors have the same color, the newly added particle 
takes this color with probability 1 — p, and the other color 
with probability p (where p is much smaller than 1). If 
the two neighbors have different colors, the new particle 
takes either color with probability 1/2. Neighbors within 
the same layer are not considered. 

Since the probability of adding a particle to a given sur- 
face site does not depend on its color, the surface grows 
exactly as with only one particle type, and is character- 
ized by the KPZ exponents \ = 1/2, and Zh = 3/2. Sim- 
ilarly, the choice of particle color at a given site is not af- 
fected by the height profile. The height profile determines 
only the moment at which a site is added, since the no- 
overhang condition requires both neighbors in the previ- 
ous layer to be occupied. If we equate layer number with 
time, domain walls move to the right or left with proba- 
bility 1/2 during one time unit, and a pair of new domain 
walls is created with probability p. This is identical to 
the Glauber model for a one-dimensional Ising chain with 



coupling J and at temperature T, with p = exp(4J//cT). 
The correlation length £ perpendicular to the growth di- 
rection is consequently £ = exp(— 2J/kT) = 1/^/p, and 
the correlation time is r = exp(— AJ/kT) = 1/p. The 
dynamical critical exponent for the order parameter is 
thus Zm = 2. Note that the "time" used for the or- 
der parameter (namely layer number) is different from 
real time, which is for each particle the moment when 
it is added to the growing surface. However, this differ- 
ence becomes negligible for sufficiently small p since the 
thickness of the surface over the correlation length, , 
is much smaller than the characteristic time, £ 2 , for or- 
der parameter fluctuations. Simulations indeed confirm 
that the order parameter and height evolve completely 
independently. A typical profile is shown in Fig. ga; the 
corresponding scaling analysis conforms to expectations, 
and is not presented here. 

The situation a > with a = b = c = can be imple- 
mented by updating sites on top of particles of different 
colors less often by a factor r < 1 compared to sites above 
particles of the same color. As the order parameter is not 
affected by the height variable, its dynamics is still the 
same as that of an Ising model, with z m = 2. The height 
profile now has domain boundaries sitting preferentially 
at its local minima, with mounds forming over domains 
(see Fig. ^o). This leads to a surface roughness exponent 
of x — 1 on length scales £, which is the case studied 
in |||. At these scales, changes in the height profile are 
slaved to domain wall motion, and the dynamic exponent 
is Zh = 2. However, on length scales much larger than 
£, the KPZ exponents of x = 1/2 and Zh = 3/2 are re- 
gained. The crossover in the roughness can be described 
by a scaling form 

([h(x,t) - h(x',t)} 2 ) = \x- x'\ 2 g(\x - x'\/0, 

with a constant g(y) for y <C 1, and g(y) ~ 1/y for y ^> 1. 

To mimic the influence of surface roughness on the or- 
der parameter (nonzero a, 6, or c in Eqs.@), the color 
of a newly added particle is made dependent not only on 
those of its two neighbors in the layer below, but also 
on the colors of its two nearest neighbors on the same 
layer, if these sites are already occupied. With probabil- 
ity 1 — p, the newly added particle takes the color of the 
majority of its 2, 3, or 4 neighbors, and with probabil- 
ity p it assumes the opposite color. If there is a tie, the 
color is chosen at random with equal probability. The 
height variable now affects the order parameter in two 
ways: (1) Domain walls are driven downhill. The reason 
is that the neighbor on the hillside of a site being up- 
dated is more likely to be occupied than the one on the 
valley side. The newly added particle is thus more likely 
to have the color on the hillside. (This corresponds to 
a > in Eq. (||).) (2) New domains are predominantly 
formed on hilltops. This is because domains on hilltops 
can expand more easily than those on slopes or in valleys, 
indicating b > in Eq. (||). Another consequence is that 
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for the same p, the correlation length £ is much larger 



than in the decoupled case, as is apparent in Figs.kLd. 




(a) 



(d) 



(b) (c) 

FIG. 2. Snapshot of the last 400 layers of simulations, for L = 200 sites, (a) The decoupled case with p = 1/90, and r = 1. 
(b) For p = 1/200, and r = 1/20, the height is coupled to the domains, but not vice versa, (c) The fully coupled case, using the 
same parameters as (b), but with updating rules that include neighbors in the same layer, (d) With r = 1, and the updating 
rules of (c), the domains are influenced by the height, but not vice versa. (Note that the profiles in (a) and (d) are identical 
since we used the same random numbers.) 
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FIG. 3. Scaling collapse of correlations and Gm in 

Fig. ^|d. For each p, the data is an average over 7500 widely 
separated layers, and for systems of size up to 8192. 

For the fully coupled case depicted in FigJ|c we find 
essentially the same scaling behavior as in Fig.^3, i.e. a 
height profile slaved to the Glauber dynamics of the do- 
mains. The most interesting case, shown in Fig.[|d, is 
when the height profile is independent of the domains 
(a = 0), evolving with KPZ dynamics, while the order 
parameter is influenced by the roughness. The dynamic 
exponent z m for the order parameter was first obtained 
by collapsing the correlation functions using Eqs. (^), as 
shown in Fig.|| These curves imply that 77 = 1, £ oc 
p-°- 542 , and r ex <f - ex 1/p, giving z m ~ 1/0.542 ~ 1.85. 



The same non-trivial value for z m is obtained by a 
completely independent measurement of the dynamics 
of domain coarsening following a quench from a "high 
temperature" (p close to 0.5) to zero temperature (p=0). 
Fig. [| shows the domain density as function of time for 
a system of size L = 16384. The resulting z m ~ 1.85, is 
in agreement with the value from the scaling collapse. 
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FIG. 4. Domain density as a function of time for 
L = 16384, averaged over 100 samples. The dotted line is 
a power-law fit (slightly shifted for better visibility) with ex- 
ponent of 1/zm = 0.542. For comparison, a power law with 
exponent —0.5 is also shown (dashed line). 

The following simple argument fails to provide the 
exponent z m c± 1.85. Consider a Langevin equation, 
x = 77 (£), for the position x of a single domain wall at 
time t. Since the motion of the domain wall is strongly 
influenced by the height profile, the noise rj(t) must have 
long-range correlations (rj(t)rj(t f )) = D\t — t f \ a , reflecting 
the dynamics of surface. This choice leads to z m = 2 for 
a > 1, and z m = 2/(2 — a) for a < 1. For a colored 
noise dominated by the slope fluctuations, a = 2/3 and 
z m = 3/2, i.e. the height imposes its characteristic time 
scale on the order parameter. This would presumably be 
the case if the domain walls were uniformly distributed 
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along the surface. However, due to their tendency to 
move downhill, they are preferentially found near valleys. 
A different scaling of the slope fluctuations in the valleys 
may be the reason for the nontrivial value of z m . Indeed, 
for short times, before the domain walls have moved to 
their preferred positions, the exponent 3/2 is seen. 

The dynamics of domain walls on a growing KPZ sur- 
face bears some resemblance to the advection of a passive 
scalar in a turbulent velocity field, which is characterized 
by nontrivial exponents and multiscaling If we ne- 
glect interactions between domain walls, and treat them 
as independent "dust particles" floating on the KPZ sur- 
face, the Langevin equation for the particle density p is 

dtp = KV 2 p + a(Vh ■ Vp + pV 2 h) + C P . (5) 

The second term describes the advection of particles 
along a velocity field v = V/i. Indeed this transformation 
maps the KPZ equation into the Burgers equation for a 
vorticity-free, compressible fluid flow ||. Equation (||) 
is a special case of Eq. (|) for m, with r = u = c = 0, 
b — a, and with a conserved noise ( p . (Together with 
Eq. (||) for the height profile, it is also a special case of 
the equations used to describe the dynamic relaxation of 
drifting polymers (9).) In the remainder, we give the re- 
sults of computer simulations for this case. The rules for 
the motion of "dust particles" are identical to those for 
domain walls. However, each particle is treated as if the 
others were not present. This means in particular that 
there is no creation or annihilation of particles. 




time 

FIG. 5. Mean square displacement of a single domain wall 
in a system of size L = 4096. The power law fit (dotted line) 
has an exponent 1.1467, corresponding to z p ~ 1.74. 

Fig. [B| shows the mean square displacement of a single 
"dust particle" in a system of size L = 4096. To obtain 
good statistics, we averaged over 512 independent and 
noninter acting particles, and used more than 40 runs. 
The best fit is obtained for z p ~ 1.74, distinct from the 
previous z m ~ 1.85, implying that the exponents depend 
on whether or not the domain walls (or "dust particles" ) 
are conserved. In contrast to the advection of a passive 
scalar in a turbulent velocity field, we find no sign of mul- 



tiscaling. Fig. ^ shows the positions of 1024 independent 
"dust particles" in a system of length L = 512. While 
there is some correlation between minima of the surface 
profile and wall positions, there are also clusters of parti- 
cles at higher elevations, indicating that particle diffusion 
is not sufficiently fast to fully adjust the density to the 
faster changing height profile. A fit of the density-density 
correlation function to (p(x)p(O)) ~ l/x 2 ^ 1 ~ Xp \ gives an 
exponent \p — 0.85. 



30 




200 400 

FIG. 6. Histogram of the positions of 1024 domain walls 

along a surface profile (indicated in grey) of size L = 512. 

In summary, the interplay between surface roughening 
and phase separation leads to a variety of novel criti- 
cal scaling behaviors. At one extreme, the height profile 
adapts to the dynamics of critical domain ordering. At 
the other, the dynamics of domain wall motion is influ- 
enced by the roughness, exhibiting new and nontrivial 
scaling behaviors. 
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